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Abstract 

Material pennittivity and permeability parameters dictate how a given material 
will react when subjected to electromagnetic fields. Simple media have been studied in 
great depth, and describe how linear, homogenous, and isotropic materials behave. 
Complex materials, those of which deviate from simple media in linearity, homogeneity, 
and isotropy, can react to electromagnetic fields in ways that differ greatly from simple 
media. 

Due to this fact, research is done to better understand these types of materials. 
The research in this study focuses on one aspect of complex media: inhomogeneity. The 
purpose is to develop an algorithm that finds the inhomogeneous permittivity profile of a 
material along its z-axis using rectangular waveguide measurements. This is 
accomplished by obtaining scattering parameters in the laboratory and using a 
Levenberg-Marquardt root search algorithm to extract the material’s inhomogeneous 
permittivity profile. 

The inhomogeneous problem is approached in two unique ways. The first method 
employs a discrete approach in which the material is assumed to be comprised of piece- 
wise constant permittivity sections. This allows the use of simpler homogeneous wave 
equations. The second method is a continuous approach which uses inhomogeneous 
wave equations to find a linear permittivity profile. 
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ELECTROMAGNETIC CHARACTERIZATION OF INHOMOGENEOUS MEDIA 


I. Introduction 

A material’s permittivity and permeability characteristics affect scattering of 
applied electromagnetic waves. These constitutive parameters influence the material’s 
behavior, whether it will act as a dielectric, magnetized material, or conductor. The 
following paragraphs provide the reader with a brief description of background concepts 
pertaining to this discussion. For a more detailed study the reader is referred to [2]. 

Simple materials are linear, homogeneous, and isotropic. These descriptions 
illustrate how the permittivity and penneability parameters behave within the material. A 
linear media is one in which the constitutive parameters are not functions of the fields, 
meaning regardless of the applied field strength, the parameters will remain unchanged. 
Homogeneity describes constitutive parameters that are constant with respect to their 
position within the material. Finally, an isotropic media is one in which the parameters 
are not dependent on field orientation. 

Complex materials on the other hand, can be nonlinear, inhomogeneous, non¬ 
isotropic, or a combination of these. Nonlinear materials have functional dependence on 
field strength, inhomogeneous profiles have varying constitutive parameters with respect 
to position, and non-isotropic characteristics define permittivities and penneabilities that 
are functions of field direction. 
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A constant dielectric profile is not always possible due to manufacturing 
inconsistencies in the material fabrication process. For this reason, a media that was 
intended to be homogeneous may instead be inhomogeneous. Averaging occurs when 
using homogeneous methods to characterize these inhomogeneous materials. This leads 
to an incorrect model when studying electromagnetic scattering and absorption. Due to 
this fact, research is necessary to enable the correct characterization of inhomogeneous 
materials. This study focuses on developing a characterization algorithm to extract the 
profile of a material that is inhomogeneous in one dimension. 

1.1 Assumptions & Scope 

Although the procedures used in this analysis could easily be extended to include 
penneability and pennittivity, this research focuses only on determining the permittivity 
profile of a non-magnetic inhomogeneous material. The study also assumes that the 
material is only inhomogeneous, and does not contain any other complex behavior. 

The problem is approached using discrete and continuous methods shown in 
Figure 1.1. In the continuous case the material with a given thickness, d is assumed to 
have a linear permittivity profile with a slope m £ , and intercept s rl . In the discrete case 
the material is described using a piecewise constant model where each section has a 
particular permittivity s n , and a total of N permittivity sections are calculated 
simultaneously. 
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Figure 1.1: Characterization of inhomogeneous material with unknown profile 

1.2 Materials & Equipment 

A number of resources were necessary to accomplish this research. The Agilent 
Technologies network analyzer was used to measure scattering parameters used in the 
profile extraction codes. Various test materials were also needed including Plexiglas, 
mica, and laminate for algorithm testing. A sample holder and a WR-90 waveguide were 
used to obtain data in the X-band frequency range of 8.2-12.4 GHz. The discrete method 
research required additional parts including a short and shim plates which were machined 
at the AFIT fabrication shop. Finally, this study required the use of facilities such as the 
microwave laboratory for taking necessary measurements. 
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1.3 Literary Foundation 

There is considerable work on constitutive parameter extraction for homogeneous 
materials. These methods typically rely on the reflection and transmission coefficients 
and offer a basis for this work. Additionally, methods have been developed to determine 
profiles of stratified media in which a single layer has unknown constitutive parameters. 

A TEM fixture was used in [9] to find the relative permittivity and permeability of 
a material. An analytical formulation of the scattering parameters was developed and 
used for this method. These scattering parameters are functionally dependent on 
frequency and time delay of the transmitted and reflected waves through a medium. This 
development provided a means for constitutive parameter extraction. 

With the advent of more powerful computers, new methods for solving these 
parameters were discovered which utilized root search algorithms [12]. This involves the 
use of scattering parameters as done before, except the material’s constitutive parameters 
are analytically difficult to solve. The root search alleviates this issue by using a 
theoretical development and experimental data to extract the desired parameters. 

The field was further developed when a direct method for parameter extraction 
was discovered [8]. This method was used to find the constitutive parameters of media 
between two known materials. This also approach was also based on a root search, and 
used an A matrix development which is closely related to the scattering parameters 
measured experimentally. This multi-layered approach provides a significant basis for 
this work, since the discrete approach in this document uses a stratified like model for 
inhomogeneous material characterization. 
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1.4 Thesis Overview 


A discrete theoretical development is formulated in Chapter 2. This theoretical 
formulation is developed for use in a discrete root search algorithm. This algorithm 
numerically extracts the discrete pennittivity profile of an inhomogeneous material using 
laboratory measurements. 

Chapter 3 describes the theoretical background for the continuous approach. 
Much like the discrete theory, this continuous theoretical framework was implemented in 
code, and a root search algorithm was used to extract the continuous permittivity profile. 

Chapter 4 contains the experimental set-up, an error analysis description, results 
found from each algorithm, and the accuracy of the two approaches. Conclusions and 
recommendations for future research are provided in Chapter 5. 


5 



II. Discrete Characterization 


The purpose of this chapter is to describe the methodology used to develop a 
discrete permittivity extraction algorithm for an inhomogeneous non-magnetic material. 
This was accomplished by modeling the material as if it were comprised of piece-wise 
constant permittivity sections. Figure 2.1 shows how this model works. If one were 
characterizing a linearly increasing profile, the discrete algorithm would produce a step¬ 
like permittivity profile description. The formulation relies on the development presented 
by Chew [4], Balanis [2], and Havrilla [7]. 

The basic theory behind the approach is first discussed using a simple two media 
interface example. The two media example is generalized to represent a stratified media 
with an arbitrary number of layers. This stratified media example is the foundation for 
the discrete inhomogeneous approach. This framework is subsequently used to calculate 
theoretical Sn parameters in the permittivity profile algorithm. 

Solving for the constitutive parameters algebraically is impossible due to the 
complexity of the problem. Therefore a root search method is used to solve for the 
desired pennittivities of each section. The root search works by comparing theoretical Sn 
parameters to experimental Sn parameters found in the laboratory using the Levenberg- 
Marquardt numerical method. This is numerical method is used calculate the unknown 
pennittivities of each section within the material, s n using 

£|[s,*0] 2 - [v;jm] 2 \<s^> n 

1 =1 
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s(z) 


E(Z) 



=0 z=d 

Discrete 


Figure 2.1: Piece-wise constant permittivity profile characterization 
In order to find each layers permittivity in this manner, a sufficient number of 
equations, L, are necessary to create a system of equations that has a greater or equal 
number of independent solutions for the number of unknowns, N. The approach used to 
obtain these equations is discussed, including the necessary hardware needed for this 
purpose. The algorithm is implemented within MatLab code to iteratively calculate the 
permittivity of each layer so that the difference between the theoretical and experimental 
values of every Su parameter is minimized to a given tolerance, 5. The Matlab default 
tolerance was used for this algorithm (10‘ 9 ). 

The discrete approach is advantageous because the underlying theory relies on a 
much simpler piecewise homogeneous treatment. However, since each layer’s 
permittivity is solved as if it were a distinct media, this requires a large number of 
independent data points. Despite this drawback, the discrete approach is significantly 
easier to implement and therefore remains one approach investigated in this study. 


7 














2.1 Discrete Theory Development 

A single-interface geometry is first discussed before investigating a multi-layered 
system. In this example two semi-infinite layers are sandwiched together and are placed 
within a WR-90 waveguide as shown in Figure 2.2. 

The first layer is free space with pennittivity and penneability of So and jio 
respectively, and extends from zero to negative infinity. The second layer has unknown 
constitutive parameters of Si and //?, and extends from zero to positive infinity. An 
electromagnetic wave is incident upon the material interface at z = 0. An artificial 
boundary denoted by a dotted line exists at z = di. This will be used to illustrate phase 
delay as the incident wave passes through the material. At the interface between the two 
materials, the energy of the wave is partially transmitted through the interface (7j), or 

partially reflected by the interface ( ). Likewise a wave originating from z=di would be 

partially reflected ( R ' 2 ) and partially transmitted ( T 2 ) through z = 0. 

The goal of the fonnulation is to represent the waves bi and c? at z=di in terms of 
the reflection and transmission of the waves bi and c; at z=0 in a way that allows material 
property extraction. After manipulation, the waves will be represented in terms of those 

found at the front and back interfaces of the material. Note that c 2 and b 2 are simply the 

waves c 2 and b 2 with a phase shift, 
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Figure 2.2 Reflection and transmission of waves through a single interface 

c' 2 =c 2 e M 

K = b ie -*^ 

K,=4ti- k l ( 2 !) 
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( 
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where, ki z is the propagation constant through material one in the z direction, k x is the 
propagation constant in the x direction, and k] is the propagation constant through the 
unknown material. Note also that a waveguide geometry is used with a width of a, and 
that the wave is a TEi 0 mode since this will be the wave excited in the laboratory. 

It can be seen from Figure 2.2 that wavehj, is the sum of the reflection of c x and 

the transmission of b 2 . It also can be seen that c' 2 is the sum of the reflection of b 2 and 
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the transmission of c x . Using these wave relationships, two equations are written that 
explicitly define the incoming and outgoing waves at z=0 


b x = R 1 c 1 + T> 2 ' 
c 2 = R' 2 b' 2 + ZJc,. 


Rewriting these equations in a more useful form an expression is obtained 


1 , R' , 

c, = —c,- -b 7 

T x T x ' 

b x =^c' 2 +(T;-^-)b[ 


( 2 . 2 ) 


(2.3) 


relating the waves on either side of the boundary. 

The reflection coefficients are related to the wave impedances of each layer 
sharing an interface. In this case layer zero has an impedance of Zn and layer one has an 
impedance of Z;. The transmission coefficient and reflection coefficients can be written 


R x = -R' 2 = 


Z\ Z 0 

z 1 + z 0 ’ 


T x =l + R x , 

T 2 = ] + R 2 . 

Using (2.3) and (2.4) the relations are manipulated to form a matrix equation 


u 

i 

' 1 

A 

c' 2 

A 

~T X 

A 

i 

A. 


(2.4) 


(2.5) 


Inserting the relations for c\ and b' 2 from (2.1) into the matrix equation in (2.5), the 
equation evolves to 


c x 

i 

~ e jk,A 

R x e~ jk ' A ’ 

C 2 


'A 

4:1 

C 2 

K 

~T X 

R x e AA 

e ~jkiA 

A. 


_Ai 

ro 

1 _ 

A _ 
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Here the A matrix is introduced. Its elements are simplified versions of the matrix 
on the left of the equation by distributing the transmission tenn to each element. The A 
matrix describes the material characteristics of layer one. This can be seen by looking 
back to the definition of the propagation constant in (2.1) and recognizing that it is related 
to the materials permittivity and permeability. Thus a formulation has been found to 
represent the waves bi and c? at z=di in terms of the reflected and transmitted waves bi 
and ci at z=0 using the properties of the material. 

Now that a simple way of characterizing a single unknown layer has been 
described, the discussion moves to a development in which there are N arbitrary layers 
within a material, each with a respective unknown permittivity and penneability. This 
development is the foundation of the discrete method. Suppose that instead of a two- 
layer system as in Figure 2.2, there were a multi-layer system as shown in Figure 2.3. 



Figure 2.3 Reflection and transmission of waves in a multi-layer system 
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From this figure a pattern quickly emerges. The problem is generalized with any 
given layer n with a thickness of d n , pennittivity and permeability of s n and u„ 
respectively, and applicable waves traveling inside and outside of its left and right 
interfaces. Also, shown in this figure is the last layer of the material ( N ), and free space 
just to the right of the material (N+ 1). 

In the previous example a single interface was investigated to develop an A 
matrix that described the waves traveling into and out of an interface containing the 
desired constitutive parameters. The multi-layer example differs from the single layer 
example in that now there are A matrices for each interface. To find the behavior of the 
entire system one simply needs to take the product of each A matrix to create an A matrix 
for the entire system. This is written as 



(2.7) 


Referring to Figure 2.3 it can be seen that the 


' i„.,— u„„ -- --n —— 


coefficient R„. The generalized value for this coefficient is 


^ Z nl 

" Z'+Z^ 


,T n =l + R„. 


This coefficient is related to the wave impedance through the layer Z„ . Also shown in 
(2.8) is the relationship between the transmission and reflection coefficients at the n ,h 
layer. 

The wave impedance can be related to the material impedance and propagation 
constants as seen in 
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(2.9) 



The material impedance and the propagation constants are then related to the 
material parameters 



V £ n 


Kn = ']°> l£ n Mn ~k x ( 2 . 10 ) 

K = 03 

The material parameters in these equations allow extraction using a numerical method. 


2.2 Constitutive Parameter Extraction 

Thus far a characteristic A matrix for each layer within the material has been 
developed. This matrix contains a description of how an incident wave will propagate 
through it. Also shown, was how to combine these A matrices to produce an A matrix for 
the entire system. These descriptions contain the permittivity, permeability and the 
thickness of each layer. From these relationships the constitutive parameters can be 
found through a numerical algorithm. 

Before this step, the A matrix for the system must be converted into an S matrix. 
This conversion relies on the relationship found in [6] 

S u S l2 

_ kJ 2l °22 

This equation allows the conversion from theoretical A parameters to theoretical S 
parameters. This step is required since S parameters are measured directly using a 
network analyzer. With the theoretical S parameters and the S parameters measured in 


1 

1 

\ 1^21 ^ 21^12 1 


1 

1 

CN 

1 


( 2 . 11 ) 
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the microwave laboratory, there now exists a means to use numerical methods for 
constitutive parameter extraction. 

In order to continue in this manner, there must be sufficient number of 
independent data points measured in the laboratory. This is accomplished by taking a 
series of measurements with a short located various distances away from the material. 

This allows the development of a set of at least N linearly independent solutions for N 
number of unknowns. 

Figure 2.4 illustrates this measurement approach within a waveguide. The 
inhomogeneous material is located at z = 0, and has a thickness of z=c!n ■ As stated earlier 
it’s assumed the material is non-magnetic for the purposes of this thesis. 

A fixed short is located a distance l\ through /„ away from the material depending 
on the distance measurement. The different distance measurements are achievable by 
inserting a waveguide shim in between the sample holder and the short. This creates an 
additional waveguide spacing between the sample and the short, providing a different 
data point. The hashed lines on the right of Figure 2.4 depict a PEC short which is offset 
from the sample holder by a given shim length l n . When calculating the system A matrix 
for each length measurement, one needs only to update the thickness for the N +1 free 
space layer by adding the applicable shim length /„. This allows for the creation of an A 
matrix for each measurement, which is used to create an S parameter equation for each 
measurement. From this point, a series of equations each corresponding to a particular 
shim length is developed to find each layers permittivity. 


14 




Figure 2.4 Sliding short measurements 
These shim lengths were calculated based on the mid-band wavelength. This 
calculation starts with the relationship between the wavelength and the propagation 
constant in the z direction 


i 2;r 
/t, — — • 
k 


( 2 . 12 ) 


The propagation constant within the waveguide in the z direction is related to the 
propagation constant in free-space and the waveguide width. 


K=Jk 0 ~ 


k 


\ U J 


(2.13) 


The propagation constant of free-space is related to the frequency and speed of light by 
the below relation 


K 


2nf 


(2.14) 
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Since the frequency at mid-band, the speed of light in free-space, and the width of 
the waveguide are known, the wavelength in the z direction can be found. To maintain 
generality in this discussion, the symbol for this wavelength is kept. 

A complete wavelength cycles through a phase of 360 degrees. It is not necessary 
to include shims through the full cycle. Any measurement beyond a half of a wavelength 
will not be independent from previous shim measurements. Therefore the shim lengths 
begin at a phase of zero and up to 180 degrees. 

The network analyzer typically needs at least five degrees of phase spacing in 
order to discern between successive measurements. So to have a confidence in the ability 
to discern independent measurements, ten degree phase spacing was used in order to 
calculate the thickness for each shim. The following formula was developed to calculate 
each shim thickness 


fllO° _ /„ 
360° X z 


(2.15) 


where / is the ri h shim thickness, n is the shim number, and X- is the mid-band 
wavelength within the waveguide. 

In this way the shim thicknesses were calculated from ten to 180 degrees. A 
“shim” of zero degrees is not necessary since it would be of zero thickness. A total of 
eighteen shims were machined. Eighteen shims plus the zero offset produced a total of 
nineteen possible measurements. 

The laboratory data obtained using these shims, and the theoretical formulation 
found earlier was coded into a root search algorithm using MatLab. The algorithm first 
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estimates a value for the permittivity of each layer (relative permeability is assumed to be 
one). These estimates are combined with each layers thickness, d,„ and used to calculate 
the theoretical A matrix for each respective layer using the theoretical relations developed 
earlier. The last layer is of free space and is produced using a respective shim offset. 

The product of these A matrices provides a representation for the entire system and 
represents a particular shim measurement. This procedure is performed in order to 
produce a system A matrix for all shim length measurements 

N +1 

<2J6) 

n =1 

The theoretical A parameters are then converted into theoretical S parameters as shown in 

( 2 . 11 ) 

DC]->[>£]• P-17) 

With the measurement data found in the laboratory, a series of equations are 
developed using the theoretical S parameters and the experimental S parameters for each 
shim. These equations are solved simultaneously to find the permittivity of each layer. 
Since the measurement construct relies on a short, only Sn parameters are measured. 

Each equation in this simultaneous computation takes the form of (2.18). Here / 
represents the equation created from the corresponding I th shim measurement. Also seen 
in this equation minimization between the theoretical Sn parameter and the experimental 
Sn parameter to a value determined by the tolerance 5 

2l[sf,(ffl,£ l ,...,£„)] ! -[s“;(®)] 2 l<<7^£,. (2.i8) 

/= 1 
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This numerical computation is accomplished using the Levenberg-Marquardt root 
search method. The details of this algorithm will not be discussed since it would detract 
from the overall discussion. However, a description of this method can be found in [3]. 
This root search method is a cross between the Newton method, and the method of 
steepest descent. First steepest descent is used to refine the estimate to obtain an 
adequate tolerance, and when sufficient accuracy is acquired the algorithm converges 
precisely to the solution using the Newton Method. 

MatLab has a built-in function that accomplishes this root search called nlinfit. A 
brief description of this function is found in the Matlab help file. Below is a line of code 
that represents how nlinfit is used 

beta = nlinfit(X,y,fun,betaO) 

where, beta is the permittivity solution to the root search, nlinfit is the Levenberg- 
Marquardt root search function, X contains the inputs necessary for the function fun that 
are particular to the problem such as frequency and material thickness, y is the solution to 
the theoretical formulation found in the laboratory data ,fun is the code that contains the 
theatrical fonnulation of Sn, and betaO is the initial constitutive parameter estimate. 

An example code using nlinfit to calculate constitutive parameters of a single 
material from two port measurements can be found in Appendix A. First, the constants of 
the geometry are assigned, these are put into the X vector. Then the laboratory data is 
read-in from a text file, this is y in the nlinfit function and represents the experimental S 
parameters. An initial guess of the pennittivity is provided to the Levenberg-Marquardt 
algorithm using betaO. Finally, the algorithm uses fun (.Appendix B) to calculate the 
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theoretical S parameters and compare those to the S parameters measured in lab to extract 
the actual permittivity (beta). 

2.3 Summary 

This chapter discussed a way to characterize layers within a material with 
theoretical A matrices. This A matrix is related to the theoretical S parameters. In turn, 
these theoretical S parameters are used with experimental S parameters to solve a series 
of simultaneous equations using the Levenberg-Marquardt algorithm. This allowed the 
extraction of the pennittivity of each layer within the material. In this way a step-like 
pennittivity profile estimate is found for the inhomogeneous material. Before proceeding 
to the results obtained using this approach, the continuous method is discussed next. 
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III. Continuous Characterization 


This chapter discusses the continuous approach for solving the inhomogeneous 
pennittivity profile. Instead of representing the material as a series of piece-wise 
constant permittivity layers, this chapter assumes the material to be a continuous linear 
inhomogeneous profile. Not all profiles of inhomogeneous materials are linear; however, 
for sufficiently thin materials a linear profile provides a good approximation. This linear 
approximation also makes it a valid approach for Radar Absorbent Material applications 
since these are typically thin. A linear profile can be seen from a Taylor series. 

, , ^e w (c), , n , . .. . s"(c)(z-c) 2 

«=o n\ 2 

If one were to approximate the function s(z) by only using the first and second tenns 

within the series. These two tenns would then represent a linear function for an 

approximate representation of e(z) at point c. 

Figure 3.1 shows the inhomogeneous material of length d 0 whose profile is solved 
in this chapter. The graph on the left of the figure shows an incident wave traveling 
through a free space medium (0), through the material with a linear pennittivity profile 
containing a z-intercept and slope (2), and finally through another region of free space 
(1). This graph describes the geometry used to calculate Su and S 2 / parameters and for 
convenience, this geometry will be referred to periodically as the forward measurement. 
This geometry is investigated first. Next S 22 and S 12 parameters are solved using the 
geometry on the right of Figure 3.1 and are associated with the reverse measurement. 
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Figure 3.1: Linear inhomogeneous profile with forward and reverse incident waves 
This graph shows an incident wave traveling through free space (1), then through the 
inhomogeneous material (2), and finally through additional free space (0). 

These theoretical S parameters are formulated in this chapter and compared to 
those measured in the laboratory for permittivity extraction. This extraction is done with 
the Levenberg-Marquardt algorithm similar to its use in the discrete case. This 
formulation closely follows the free space development found in [11], except that a wave 
is propagating within a WR-90 waveguide environment. 

3.1 Forward Measurement Field Descriptions 

This discussion assumes an incident wave originating from the -z direction. The 
forward transmission and reflection coefficients will contain a subscript/to distinguish 
them from the reverse coefficients developed later in this chapter. 

The electric and magnetic fields must be found before constructing theoretical 
transmission and reflection coefficients. First a formulation for the electric field within 
region zero is developed. The y component of the electric field is a function of the 
position x and z. The field’s dependence on z is denoted by the psi function as shown 
below 

E y {x,z) = sm(k x x)y/{z), (3.1) 
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where, k x = —. 

a 

The electric field in region zero is written directly by inspection. From the figure, 
forward and reverse traveling waves are expected within the region since there are both 
incident and reflected waves. For this reason, 

y/, = R f e ik - z +e jk ' z , (3.2) 

where, ^represents the electric field’s z dependence for region zero, and k z is the 
propagation constant in the z direction. This propagation constant is defined as 

K = -\kl-kl 
K ~ o ’ 

where, k x is the x-di reeled wavenumber, ko is the free space wavenumber, and a denotes 
waveguide width. 

There is only a transmitted wave \j/ x in region one. So this is also written by 
inspection 

y/ x = T f e jkM( '~ z) . (3.3) 

Now that the electric fields for region zero and one are solved, the magnetic fields 
are found easily from Faraday’s Law 

H x0 (x,z) = sm(k x x) 7 d ^° . (3.4) 

(Oju tj dz 

Although there is also a z component of the magnetic field, it is not needed for 
enforcement of boundary conditions, and therefore only the x component of the magnetic 
is used. The magnetic field for region zero simplifies to 
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H t0 (x, z) = sin(M) k A ' ~ ] 


(3.5) 


Using this same procedure the x component of the magnetic field for region one is 

H xi = (x, z) = sin (k x x)-^—T f e jkz(da ~ :) . (3.6) 

co/u 0 

The process to find the electric and magnetic fields for region two is more 
involved since the region is inhomogeneous. A differential equation is derived from 
Maxwell’s equations in order to describe the electric fields in region two. 

Faraday’s law shows that the electric field will circulate around a magnetic field, 
and Ampere’s Law shows that a magnetic field will circulate around an electric field 

V x E = -jcojuH (3.7) 

V x H = jcosE. 


In these equations E represents the electric field, H the magnetic field, e and // the 
pennittivity and permeability, and co represents the angular frequency of the wave. 

Expanding the curl operator and recognizing y invariance since the material 
properties do not change in the y direction, Faraday’s law yields 


-x —-+ v 
dz 


^ =-ja J u{xH x + yH v+ zH z ). 


By equating the vector components a set of three equations results 


-r-- = j(OMH x 
dz 



T 1 = j^H v 
dz 


dE 

— =-jco/iH z . 

dx 


E y ,H x : TE-- 
E z ,E x ,H x : TM z 
E y ,H z : TE-’ 


(3.9) 
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Ampere’s law results in a similar fonnulation 


-x—-+ v 
dz 


dz dx 


= jos(xE x + yE v + zE z ). 


(3-10) 


Again, a set of three equations is found by equating the vector components from (3.10) 


dH v 

— = -j(osE x 
dz 



= -jcosE 
dz 


dH 

—- = jcosE z . 

dx 


E x ,H y :TM Z 
E y ,H x ,H r :T\E 
E_,H V : TM r 

z 


(3.H) 


These results make physical sense because in the TE Z case the magnetic field 
circulates around the electric field, and in the TM Z case the electric field circulates around 
the magnetic field. 

The problem geometry dictates y invariance. This means that the equations are 
decoupled, and a TE Z excitation will not couple into a TM Z field set. This simplifies the 
problem considerably because the TE Z field excited in lab will not generate TM Z fields. 

As a result only TE Z waves are required in this analysis. 

In particular the problem will consist of the TE[ 0 mode, since this specific mode 


will be used within the waveguide. Figure 3.2 shows the expected waves within the 
guide. 


24 



te:, 


TEE, 


te:. 


te:. 


te:. 


Figure 3.2: TEio modes in WR-90 waveguide 


Solving the TE Z equations from (3.9) in terms of H x and //-, and inserting them 
into the second equation in (3.11), a partial differential equation in terms of the electric 
field is derived 

d 2 E d 2 E 

—^ + —^ + k 2 (z)E y = 0. (3.12) 

ox oz 

This is a homogeneous partial differential equation that represents the electric 
field in the inhomogeneous material. Note that electric field E y is a function of x and z 
for region two and is written 

E v (x, z ) = sin (k x x)y/ 2 (z). (3.13) 

By assuming separation of variables (3.12) is reduced to simple derivatives 

-k 2 x E v + ^f + k 2 (z)E y =0. (3.14) 

dz 

Note that the partial double derivative of E y with respect to x is just -k x 2 E y since there are 
only standing waves in the waveguide along the x-axis. 
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This result simplifies further using 


k 2 = k] + k; (3.15) 

to obtain 

d - E f + k 2 (z)E v = 0. (3.16) 

dz 

Where, k z is functionally dependent on position z because the material is inhomogeneous 
in this direction, and k v is zero since the material is y invariant. 

There is no general analytic solution to this differential equation so a particular 
profile is assumed. For sufficiently thin materials a linear profile can be assumed without 
losing significant accuracy. This assumption dictates a linear propagation constant 

k: (z) = k 2 (z) - k; = co 2 £ 0 n 0 s r (z) - k; (3.17) 


This simplifies further into a more useful representation 


k 2 s r (z) - k; = k 2 (s rl + m £ z - -f). 


(3.18) 


From (3.18) the pennittivity profile is much easier to discern. By factoring out 
the propagation constant of free space the equation is left with the relative pennittivity 
intercept £ r j and the relative pennittivity slope m £ of the inhomogeneous profile. If the 
slope were to go to zero, the resulting equation would be exactly equal to that for a 
homogeneous material. Substituting (3.18) into the propagation term in (3.16) yields 

■ E + klle r , + m,z -1i]L,(r) = 0 (3.19) 

[dz~ k 0 J 


where, the notation used for the electric field in region two is fS. 
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The differential equation in (3.19) is in a complicated form. A change of 
variables is used to simplify the problem 


77 = — \m e z + i m-tt , 


(3.20) 


or equivalently, 


1 ( m 


tj- \ s 


(3.21) 


One could substitute the double derivative with respect to z with the double 
derivative with respect to 77 by knowing 


dy/ 2 {z) _ dy/ 2 ( 77 ) dr/ 
dz drj dz 


(3.22) 


d 2 y/ 2 (z) _ d f dy/ 2 (z)~\dri 


dz drj\_ dz \ dz 


(3.23) 


The derivative of 77 with respect to z is 


dd 1/3 j 2/3 


1/3 7 2/3 

= "T K » 


(3.24) 


so that (3.19) becomes 


-A«C) 2 + i 0 ! K, +m ,U. A 

clTi m s ^ /c 0 


f k ; 

n- £ r\ ~~Tl 


-Jj] ¥1 ( 7 ) = 0 - 

/Cr, 


(3.25) 


Simplifying this expression results in the differential equation 


t + ?7 ^( 77 ) = 0. 


(3.26) 
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The solution to this differential equation is found using the Frobenius method. 

One will find the solution to be the scaled sum of two fractional-order Bessel functions as 
seen in [10] 

1 ^2O) = w 1/3 [AJ 1 /3 (w) + BY 1/3 (w)], (3.27) 


where, A and B are constants, and an additional change of variables was made 



2 K l 

w = - — (m e z + £ r 

3 m. 



Bessel functions make physical sense, since the waves traveling forward and 
backward should appear sinusoidal. Also, one would expect the waves to diminish in 
amplitude and decrease in wavelength as they travel through the material to account for 
loss and changing permittivity within the material. This phenomenon can be seen from 
the relationship 


A = — = —^=, 

f 

where v and/are the velocity and frequency of the wave respectively, and c is the speed 
of light in free space. As pennittivity increases the wavelength would naturally decrease. 
Since these relationships intuitively represent what is expected in the problem, the 
solution in (3.27) is reasonable and the formulation continuous to finding the magnetic 
field in region two. 
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The magnetic field can be found as previously accomplished using Faraday’s 
Law. In this case, a change of variables was introduced, and therefore the derivative of w 
with respect to z must be added to the equation 

H x2 (x,z) = sin(k x x) J d V 2 dw (3.2 8 ) 

a>jU 0 dw dz 

Using the property of Bessel functions below, taking the derivative of the electric field 
becomes quite simple 

^[ zV &( z )] = zV &-i( z )- 

This holds true for any Bessel function % v (z ). Using this equation the derivative of the 
electric field is 


^- = w m lAJ^ a (w) + BY_ m (w)]. 

dw 


(3.29) 


The derivative of w with respect to z is 


— -k 
, K 0 

dz 


m e Z + £ rl-Jl) 

V K J 


1/2 


(3.30) 


Using (3.29) and (3.30) the magnetic field for region two becomes 


H x2 (x, z ) = sin(k v x) —— k 0 w 1 ' 3 


m £ z + s rl - 

V J 


[AJ_ m {w) + BY_ m {w)\. (3.31) 


Now that both the electric and magnetic fields in all regions are known, boundary 
conditions are applied to find the forward measurement’s reflection and transmission 
coefficients. 


29 



3.2 Forward Measurement Transmission & Reflection Coefficients 

Since the tangential electric field components are continuous across the boundary 
at z=0, the electric field in region zero and region two must be equal 

^o(0) = ^ 2 (0). (3.32) 


Note that the electric fields x dependence cancels and only the z dependent term 
remains. This cancelation will occur when applying all boundary conditions including 
those for the magnetic fields. Using the equations for the electric fields in regions zero 
and two, and evaluating at z=0, results in 

Uo(°) = R f + 1 

y 2 (0) = wf [AJ m (w„)+BY v ,(w a )\. 


Equating these two fields at the boundary yields 

1 + R f = W T [AJ+ BY m (w 0 )\, 

where, 


(3.34) 


f , V 

\ m sJ\ 


K j 


3/2 


The next boundary is z = do in which the electric fields are equal in regions one 


and two 


Wx (d G ) = y 2 (d a ). (3.35) 

At this boundary the electric fields for both regions are 

Vi (d 0 ) = T f 

(3.36) 

¥i(d () ) = w, 1 ' 3 [AJ m (w x ) + BY m (w,j], 
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where, 


2 f k n Y k 2 

W l = T - ™A +£ r l-TJ 

v y V 


Equating these tangential electric field values according to (3.35) yields 

r / =w 1 1/3 [yU 1/ 3(w 1 ) + 5^ /3 (w 1 )]. 


(3.37) 


Next the boundary conditions for the magnetic fields are applied, starting with the 
continuity of tangential magnetic fields across regions zero and two using z = 0 


H xO (0) = H x2 (0). 


(3.38) 


Equating the two formulations at this boundary the following equation results 


{ - Rf= lJLJL_ [AJ^ m (w 0 ) + BY_ m (w 0 j\. 


(3.39) 


Finally, the continuity of the magnetic field across regions one and two at z = do is 
applied 


H x M = H x2 (d»). 


(3.40) 


Using this boundary condition and the equation for the magnetic field in region two 
developed earlier (3.30) the last of the four boundary value equations for the forward case 
is found 


m A + e r i-tt 


[AJ_ m (w,) + BY_ m (w,)\. 


(3.41) 


At this point four equations have been developed from the boundary conditions. 


These equations define the reflection and transmission coefficients for the forward case. 
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Note that the constants A and B must be removed by solving these equations 
simultaneously 


where, 


1 + Rf ~ W T [^l/3( W o) + ^^/3( >V o)]’ 
T f =w^[AJ m {w l ) + BY m (w l )\, 

1 ~R f = jC\ f w o 3 [^- 2 / 3 K) + BY_ m {w 0 )\, 
T f = 7C 2/ < [AJ_ m {w x ) + BY_ m {w x )], 


(3.42) 


r 

1 / 

r 

2 / 


= K 
K 

= K 
K 


k 1 

( £ 

V^rl ,2/ ’ 

K 0 

(m £ d 0 + s rl - 



Transmission and reflection coefficient equations are found by solving the set of 
relations in (3.42). The coefficients are again labeled with a subscript/to emphasize that 
they are for the forward measurement; those taken by an incident wave originating from 
the —z direction 


|V 

1/3 

N l f + N 2f 



D\fD 2 f + ^3 f^4f 


(3.43) 


R f = 2 


N if + N 4f 


D\/D 2 f + 


- 1 , 


(3.44) 
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where, 


N i f = J m (' w i) [ Y m (' w i) - jC 2 ;/ A/ 3 Ofi )] 

A/ = A(ffi)[./A/A/3(*fiWl/3( W l)] 

A/ = ^l/3( W o)[-^/3( W l)~ i^2/-^2/3( W l)] 

A/ = 7 l/3 Of) ) [i A /A/3 Al )~J 1/3 Al)] 

A/ = \jm A) - J C 2 f Y -in A )] 

A/ = [y Ay J~m Ao) + *A/3 Ao)] 

A/ = [A Ao) + yA/A/3 Ao)] 

A/ = [y'A/A 2 /3(wj)— 

Fortunately the solution for the transmission coefficient can be reduced using the 
relationship for any Bessel function E , v , 

- 7 - [4v (« z )] = -«& + i (« z ) + - 4 (« z )> (3-45) 

az z 

and the identity, 

J v ( z )Y v '(z )-Y v (z )J[ (z ) = —, (3.46) 

nz 

which is also kn own as the Wronskian. This allows the transmission coefficient to 
reduce 

( -4 jC 2f V -i-i 

T ,= A L [A/Az+A/A,] ■ (3.47) 

AAl>V 

There were no such simplifications that could be made for the denominator term or for 
the numerator for the reflection coefficient and therefore the form shown in (3.44) 
remains. 
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The transmission and reflection coefficients are related to S parameters as shown 

below 

Sff(a),£ rl ,m e ) = R f (G),e rl ,m e ) 

.4o ) 

S 21 (a>,£ rl ,m s ) = T f (co, e rl ,m £ ). 

Next the theoretical S parameters for the reverse measurement are developed. 

3.3 Reverse Measurement Field Descriptions 

The development in this section closely follows the approach used in section 3.1. 
In this section however, the wave originates in the z direction, as shown in Figure 3.3. 

The reflection and transmission coefficients belong to the reverse measurement set and 
are denoted with a subscript r. These coefficients are directly related to the S 22 and S 12 
scattering parameters. 

The formulation begins with finding the electric field within region zero. Using 
the same convention for the electric field, i// 0 can be written by inspection as 

V* =T r e jk -\ (3.49) 

A reflected wave and an incident wave will be present in region one so that y/\ is 

y/ x - + Rr6 jkM 0 -z)_ (3.50) 

Now that the electric fields are known for regions zero and one, the magnetic 
fields are found from Faraday’s law 

H x0 (x,z) = sm(k x x) J d ^° . (3.51) 

coju tl dz 

Therefore the magnetic field for region zero is 
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(3.52) 


H x0 (x,z) = sin(A: v x) —— T r e jKz . 

O)/J 0 

In a similar fashion the magnetic field for region one is 

H xl (x,z) = sin(k x x)^— \ R r e ik * (d °~ z) - e~ jk * {d °~ z) ~\. (3.53) 

' (y// 0 L J 

At this point the electric fields and magnetic fields have been found for region 
zero and one. For region two, no additional formulation is needed since the wave within 
this region will behave similarly whether the incident wave originates from the -z or z 
direction. For this reason, (3.27) and (3.31) can be also used in the reverse development. 


3.4 Reverse Measurement Transmission & Reflection Coefficients 

The tangential electric field is continuous across the boundary z=0. Because of 
this, the electric field in region zero is equal to the electric field in region two as noted in 
(3.32). At this boundary the z dependence of the electric fields for region zero and region 
two are 


wM = T r 

¥ 2 (0) = wf [AJ m {w 0 ) + BY m {w 0 )\. 


(3.54) 


Note that wo is equivalent to that introduced for the forward case. Equating these two 
electric fields results in 


T r = w 1 ' 3 [AJ m (w Q ) + BY m (w 0 ) . 


(3.55) 


Next, the boundary condition for z=do is applied. Here the tangential electric 
fields are equal for regions one and two. The fields for the two regions are shown below 

¥\ (d 0 ) = \ + R r 

¥ 2 ( d o) = W T [A/^Nql + ZN^lvv,)]. 
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Note again that w/ was defined earlier in (3.36). Equating these two electric fields results 
in the relation 


1 + R r =w\ ,2 [AJ m (w l ) + BY m (w l )\. 


(3.57) 


Next boundary conditions are used to find two additional equations from the 
continuity of tangential magnetic fields. This development starts with the magnetic fields 
across regions zero and two at z = 0. At this boundary the magnetic fields are 

-kX 

(3.58) 


H x0 (x, 0) = sin(k x x)- 


H x2 (x,0) = sin(k x x) 


0)/J Q 

jk,Xo 


. 1/2 


( 'Xu 


■ r \ ,2 

V K 0 J 


2/3 Oo)]- 


Equating these two magnetic fields across the boundary a third equation results 

k 2 . 


T _ ~jk 0 w Q 


1/3 f i 2 V /2 


£ " 1 P 
V *o J 


[Aj -2/3 (w*)) BY 2/i (w Q ) • 


(3.59) 


Finally, the magnetic fields across the boundary z=do are also equal. These 
magnetic fields are 


kf xt (x,d () ) = sin (k x x)—*-(R, - 1) 

m, 


H x2 (x,d 0 ) = sin(^x) 


jK w \ 


1/3 f k 


(DjU 0 


m e d 0 + £ r i-yy 

V J 


1/2 


(3.60) 


[AJ_ m (w,) + BY_ m (w,)\. 


This results in the last equation 


R -1 = 


A w i 


1/3 


f 




, p 

-TT 
^0 J 


1/2 


\aj i )+BY 2/jiw^. 


(3.61) 


The four equations describing the transmission and reflection coefficients for the 
reverse measurement are written as 
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(3.62) 


T r = w o 3 [ AJ 1/3K) + BY m Ao )] 

1 + K = [ AJ mi w \) + BY \^ w \)\ 

T r = jC h .wf [AJ 2/3 (h> 0 ) + BY_ m {w 0 j\ 
R r ~ l = jC 2r w\ /3 [A/_ 2/3 (h',) + flK 2/3 (W|)] 


where, 


C„ 


a 



k 

= ~f~(m e d 0 +s rl 

K 



With some algebra the constants A and B drop out and relations for the 
transmission and reflection coefficients are 


U1 

1/3 

A+A,- 



_ArAr + ArAr . 


(3.63) 


where, 

N u- = J m 0 % ) ft/s ('% ) - J C \J-m (% )] 

Ar = A ( W 0 ) [i A A-2/3 ( W 0 ) _ A/3 ( W 0 )] 
Ar = ^U/3 C W 1) [y 2/3 C) *^1/3 C W 0)] 

A r =[YvM)-jC 2 J^{w x )\ 

A r = [iAr^A^Ao) ~ A/3A0)] 

A. ^AAoWAA^A,)] 

A r = |A/3 A| ) — 7 A A-2/3 Al )] 


(3.64) 


As done for the forward case, the transmission coefficient can be further simplified using 
(3.45) and (3.46) to yield 
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T = 


' -4 jC u 4 


2 . \l/3 


n{w Q w x ) 


[AA+AAf 1 ’ 


(3.65) 


Now that the reflection and transmission coefficients are found for the reverse 


measurement, the associated theoretical S parameters are simply 


Sn (®> > m s ) = £ n » ) 

(co, £ rX ,m £ ) = T r (co, £ rX ,m £ ). 


(3.66) 


3.5 Numerical Computations of Slope and Intercept 

To summarize, the four transmission and reflection coefficients are related 
directly to S' parameters as shown below 

T f (co, £ rl ,m £ ) = S* (co, £ rl , m £ ) 

R f (co, £ rl , m £ ) = S'* (co, £ rl , m £ ) 

T r (co, £ rX , m £ ) = S* (co, £ rX , m e ) 

R r (co, £ rX ,m £ ) = S* (co, £ rl ,mj. 

Before continuing, a quick check is done to see if the mathematical representation 
of scattering parameters reflects physically expected results. This allows a check on the 
accuracy of the formulations. Unfortunately, allowing the slope m e to go to zero to verify 
the solution reduces to a homogeneous relationship does not result in a useful test. This 
is because as m s approaches zero, fields would become infinite within the inhomogeneous 
region. Even taking the Bessel functions to their asymptotic approximations did not 
allow for a testable result. 

However, some other checks can be made. Since the inhomogeneous profile is 
asymmetric, the equations for the forward and reverse reflection coefficients should not 
be equal. Referring back to (3.44) and (3.46) it is apparent that these coefficients indeed 
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differ. The theorem applying to linear and isotropic materials can be used. For both 
homogeneous and inhomogeneous materials reciprocity signifies that the forward and 
reverse transmission coefficients should be equal. Referring back to (3.47) and (3.65) it 
is seen after inserting the values for C//, C^/, C;,., C>, wy, and vv/, that these transmission 
coefficients are indeed equal. This provides reasonable confidence that these scattering 
parameter fonnulations are accurate. 

The four scattering parameter equations are solved using the Levenberg- 
Marquardt algorithm to find the relative permittivity z-intercept and slope within the 
inhomogeneous material 

El[«,*’(»■*„■'».)]' -[S“ p (®)] 2 1 <S^e rl ,m„ (3.67) 

i 

where the scattering parameters are indicted using 7=21,11,12, and 22. 

3.6 Summary 

In this chapter theoretical scattering parameters were solved allowing a linear 
pennittivity profile extraction. This was accomplished by splitting the geometry up into 
three regions: zero and one of homogeneous free-space, and the second of an 
inhomogeneous material. From this point, boundary conditions were applied to solve for 
the forward and reverse incident wave’s reflection and transmission coefficients. These 
coefficients were directly related to the scattering parameters. 

Finally, as done in the discrete case, an algorithm was developed via MatLab to 
numerically solve for the linear profile using the Levenberg-Marquardt algorithm with 
experimentally captured scattering parameters. 
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IV. Analysis & Results 


This chapter discusses the steps taken in the laboratory to capture data, and the 
results from the discrete and continuous pennittivity profile algorithms. A short error 
analysis section shows how pennittivity profile uncertainty was calculated. All results 
are presented using Matlab plots to clearly show the accuracy of each algorithm in 
various testing scenarios. Also included in this chapter is a comparison between both 
approaches and a discussion of advantages and disadvantages of each. 

A description of the equipment used, the materials tested, and the measurements 
taken are included for completeness of the discussion. The discrete algorithm required a 
significant number of measurements, and so details on data gathering are included for this 
approach. The continuous approach required only one measurement, and so only a brief 
description is required. 

This chapter begins with the experimental set up, the calibration and measurement 
procedures, error sources, algorithm results, and concludes with a summary. 

4,1 Experimental Set Up 

The scattering parameters were measured using the Agilent Technologies 
Network Analyzer shown in Figure 4.1. Shown is the S Parameter measurement screen, 
the two ports that connect to the waveguide, and the keypad. 

Several waveguide parts were needed to collect the necessary data. Shown in 
Figure 4.2 are the short (A), the sample holder (B), and material samples, including 
Plexiglas (C) and Mica (D). 
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Figure 4.3 shows a picture of the shims displaying their thicknesses, which were 
calculated earlier. Starting from top to bottom, and going from left to right, the shims get 
progressively smaller. The last shim shows the face of each offset. In table 4.1 the 
specific thicknesses of each shim is listed. 

Finally, Figure 4.4 and Figure 4.5 show the experimental set up for the discrete 
and continuous measurements, respectively. Figure 4.4 shows the coaxial cable (A) that 
connects the network analyzer to the waveguide (B), the sample holder (C), the shim used 
for the applicable measurement (D), and the short (E). Note also in this figure that only 
measurements were taken in the discrete case using one port (Sn). Figure 4.5 shows a 
slightly different set up. In this figure measurements were taken using two ports (Sn, S 21 , 
S 12 , S 22 ), both labeled A, where the left was port one, and the right was port two. The 
waveguide is also shown (B) along with the sample holder (C). 

Table 4.1: Shim thickness measurements 


Shim 

A 

B 

C 

D 

E 

F 

G 

H 

I 

Thickness (mm) 

19.87 

18.76 

17.60 

16.54 

15.43 

14.34 

13.23 

12.10 

11.02 

Shim 

J 

K 

L 

M 

N 

0 

P 

Q 

R 

Thickness (mm) 

9.96 

8.79 

7.71 

6.61 

5.52 

4.40 

3.31 

2.19 

1.10 
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Figure 4.1: Network analyzer 



Figure 4.2: Short, sample holder, and materials 
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Figure 4.3: Shim offsets 



Figure 4.4: Single port measurements for the discrete approach 
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Figure 4.5: Two port measurement for the continuous approach 


4.2 Calibration and Measurements 

For both discrete and continuous measurement data sets, TRL calibration was 
used to ensure data accuracy. A brief description of TRL calibration theory is found in 
[6]. The network analyzer’s built in TRL calibration kit was used prior to data collection. 

A total of thirty eight measurements were taken for the discrete case when testing 
on a stratified material. This was done to find additional data points, and can only be 
accomplished for media that are sufficiently thinner than the sample holder to make these 
measurements possible. The stratified configuration is the only material that was tested 
in this way. 

The first nineteen measurements were taken with the sample to the left of the 
sample holder. The next nineteen measurements were taken with the sample to the right 
of the sample holder. Figure 4.6 shows the two different measurement configurations. 
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The left side of the figure shows the sample to the left of the sample holder, and the right 
side of the figure shows the sample to the right of the sample holder. The figure also 
shows the waveguide (B), sample holder (C), shim (D), and short (E). In each set of 
measurements the short immediately followed the sample holder. 

The first measurement was taken with no shim, and then measurements continued 
with the smallest shim (R in Figure 4.3). Each subsequent measurement was taken with 
an increasingly larger shim. The last measurement being the final shim (A in Figure 4.3). 
The first test of the discrete algorithm was using a stratified media with these thirty-eight 
measurements. This stratified media was comprised of a thin sample of Plexiglas, 
laminate material, and mica. At the end of taking all discrete measurements there were a 
total of thirty-eight data files, each containing the Sn scattering parameters associated 
with a given frequency. For the pseudo continuous inhomogeneous material, only 
nineteen measurements were possible since the material was too thick to create a right 
justified measurement set. 

For the continuous case one two-port measurement was needed. In that instance 
the shims were not used. For these measurements an inhomogeneous material was 
constructed to collect the data as well as a homogeneous material for initial testing. In 
this instance, one file was created with Sn, S 2 i, S 12 , and S 22 scattering parameters 
associated with a given frequency. 
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Figure 4.6: Left justified and right justified sample measurements 


4.3 Error Analysis 

The error analysis used in this section is similar the differential treatment in [5]. 
There are numerous sources of error when using these algorithms because they rely on 
experimental data for conducting the root search. This can be seen by revisiting the 
extraction formulae for the continuous case 

Xl [st y (co,s rX ,m £ )^ -[^(ru)] 2 1< S -A s rl ,m E , (4.1) 

i 

and the discrete case 

i=i 

Not all errors are included in this discussion since many of them are difficult to 
measure and present minimal uncertainties. Examples of such errors include excessive 
movement of the coaxial cables, misalignment of the measuring constructs shown in 
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Figures (4.4) and (4.5), and uncertainties internal to the network analyzer. These are not 
investigated in this analysis. 

Other errors can be approximated more easily and contribute more to the overall 
results and are used in the analysis. These include inaccuracies in material thickness 
measurements, placement of the material within the sample holder, and the uncertainty in 
the shim lengths for the discrete algorithm. The difference between the discrete and 
continuous case is that the continuous algorithm doesn’t rely on shims and therefore has 
fewer sources of error in comparison. 

Since the discrete case presents more error possibilities it is described in detail. 
The continuous case would follow a similar treatment without the shim measurement 
uncertainty. 

The complex permittivity profile of the material can be written functionally as 

s = s’{m, p, 1 ) - js"(m, p, /). (4.3) 

This permittivity is comprised of a real component s' and an imaginary component s”. 
These are both functions of the material thickness m, the position of the material within 
the sample holder p, and the length of each shim /. There is uncertainty in each of these 
measurements leading to a perturbation from the actual permittivity profile 

A s = s'(m 0 + Am, p 0 + Ap, l 0 + A/) - js"(m 0 + Am, p 0 + Ap, 1 0 + A/). (4.4) 

This equation shows that the uncertainty of the pennittivity is a function of the 
actual measurements for the material, placement, and shim length denoted with the 
subscript 0, and the uncertainty of these same measurements shown by Am , Ap, and A/. 
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The real permittivity could be written using a Taylor series as 


s'(m 0 +Am,p 0 + Ap,l 0 +Al) = s\m 0 ,p 0 ,l 0 ) + 

d£'(m 0 ,p 0 ,I 0 ) Am | de'(m 0 ,p 0 ,l Q ) | de'(m 0 ,p 0 ,l 0 ) A/ | (4.5) 

dm dp dl 


This is also true for the imaginary permittivity; however, to simplify the discussion only 
the real pennittivity will be discussed. The imaginary permittivity uncertainty follows a 
similar process. 

If the error is small, the higher order terms can be ignored in the Taylor series 


s\m 0 + Am , p 0 + A p, / 0 + A/) - s'(m 0 ,p (j ,l 0 ) ~ 

de'(m 0 ,p 0 ,I 0 ) Ain | ds'(m 0 ,p 0 ,l 0 ) | de'(m 0 ,p 0 ,l 0 ) A/ (4.6) 

dm dp dl 

The partial derivatives can also be approximated if the errors are sufficiently small 


_ gKdU) Am x e'(m 0 + Am,p 0 ,l 0 )-e'(m 0 ,p 0 ,l 0 ) ^ 
dm Am 

_ d£'(m 0 ,p 0 ,l 0 ) g e'(m 0 , p 0 + Ap,l 0 )~ £'(m 0 , p 0 ,l 0 ) 
p dp Ap 

A „> _ d£\m 0 , p 0 ,4) A/ ^ e'(m 0 ,p 0 ,l 0 + A/)-g'(ffl 0 iPM_ m 
7 5/ A/ 


(4.7) 

(4.8) 

(4.9) 


The magnitude of each uncertainty is found and combined to find the total uncertainty in 
the resulting permittivity A s' extracted by the algorithm 

| As' |=| s'(m 0 + Am, p 0 + Ap, l 0 + Al) - s'(m 0 ,p 0 ,l 0 ) H A s' m \ + \ As' p \ + \ As\ \. (4.10) 

The absolute values are used so that the errors associated with each uncertainty do not 
cancel with each other. This results in a total permittivity error that presents a worst case 
scenario regarding measurement errors. 
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The uncertainties Am , A p , and A/ can be positive or negative. For small 
uncertainties the resulting permittivity perturbations will be symmetric about 
e'(m 0 ,p 0 ,l 0 ). As the uncertainties increase the positive and negative errors can become 

significantly different from each other, and need to be calculated separately. It is 
assumed in this analysis that this difference is sufficiently small and does not warrant 
separate calculations. In this manner the uncertainties were calculated for the real and 
imaginary permittivity profiles for the discrete and continuous cases. 

In the discrete case A/ was calculated by finding the standard deviation from the 
shim measurements, and using this uncertainty per shim during permittivity extraction. 
The material measurement errors (Am ), for the discrete and continuous cases were 
calculated using the standard deviation among material measurements, and depended on 
what material was being tested. The material placement error within the sample holder ( 
A p ) was set to two tenths of mm for the calculations. These uncertainties will be seen on 
the graphs with the use of error bars. 

4.4 Algorithm Results 

This discussion starts with the results from the discrete case. These first graphs 
compare a material as measured within the stratified configuration (thin Plexiglas, 
laminate, and mica) with the inhomogeneous algorithm (red) and compare them with the 
same material that is measured individually using a homogeneous algorithm (blue) with a 
two-port scheme. Also note that each graph contains error bars for both the 
inhomogeneous and homogeneous algorithms using their respective color. 
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In Figure 4.7, Figure 4.8, and Figure 4.9 we see two graphs per figure. The first 
plots the real part of the relative permittivity versus frequency for both the 
inhomogeneous and homogeneous algorithms. The second plots the imaginary relative 
permittivity versus frequency for both the inhomogeneous and homogeneous algorithms. 
The figures are displayed in the order in which the materials are configured within the 
sample holder in the stratified configuration. 

The first material encountered in the stratified media is Plexiglas. The real 
relative permittivity for Plexiglas is close to the actual value, but does differ at higher 
frequencies. The imaginary relative permittivity produced by the algorithm is strongly 
consistent with the actual imaginary relative permittivity. Immediately one can see the 
result error can play in the pennittivity extraction using the discrete algorithm. 

The second material in the stratified media was laminate. In the corresponding 
figure, a strong agreement with the imaginary results is seen, however again a slight 
divergence between the real permittivity in the inhomogeneous algorithm and the real 
permittivity in the homogeneous algorithm is seen. 


Plexiglas: Real Relative Permittivity vs. Frequency 


Plexiglas: Imaginary Relative Pennittivity vs. Frequency 



Figure 4.7: Real and imaginary results for Plexiglas 
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Laminate: Real Relative Permittivity vs. Frequency 


Laminate: Imaginary Relative Permittivity vs. Frequency 
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Figure 4.8: Real and imaginary results for laminate 


Mica: Real Relative Permittivity vs. Frequency 
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Mica: Imaginary Relative Permittivity vs. Frequency 
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Figure 4.9: Real and imaginary results for mica 
Finally, the third material encountered in the stratified sample is mica. There is 
again a slight divergence in the real permittivity results, and strong correlation in 
imaginary permittivity results. 

Since there are over thirty-eight measurements in the discrete case, the cumulative 
effect of all uncertainties begins to increase the amount of error possible in the results. 
This can be easily seen when comparing the error bars between the inhomogeneous and 
homogeneous results. For this reason, great care must be taken when measuring the 
scattering parameters for the discrete case. This complex measurement environment is 
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believed to be the cause of the inaccuracies in the discrete algorithm. The smaller error 
seen in the imaginary calculations is believed to be caused by the fact that the material 
has relatively low loss. If the material has greater loss, it is likely the error bars would be 
just as large in the imaginary calculations as they are in the real calculations. 

Next the results for the continuous case are discussed. This algorithm was tested 
with a homogeneous Plexiglas material and later with an inhomogeneous material created 
using layers of blue heavy-weight paper and yellow plastic. 

The first test was a Plexiglas material to get a sense of accuracy with the 
algorithm. Since the material in region zero has a relative permittivity of one (free space) 
the z-intercept was assumed one, and only the slope was computed. For this comparison, 
the inhomogeneous slope was calculated with the continuous inhomogeneous algorithm 
and was plotted against the rise over run permittivity of the Plexiglas (calculated using 
homogeneous code). Figure 4.10 shows these results. Note that the slope is represented 
in these graphs as relative pennittivity divided by distance in mm. 

It can be seen from this figure that the real relative slope is fairly accurate, with a 
slight divergence above 11.5 GHz, and for the imaginary relative slope there is a smaller 
divergence. Since the mathematical solution of the differential equation in region two 
has an instability when slope is zero, a completely homogeneous material could not be 
used to test the algorithm to find both intercept and slope. Instead an inhomogeneous 
material was made for further testing. 


52 



Plexiglas: Real Relative Slope vs. Frequency 



Frequency (GHz) 



Frequency (GHz) 


Figure 4.10: Real and imaginary relative slopes for Plexiglas 
Next the continuous algorithm was tested with a pseudo continuous permittivity 
profile by interspersing layers of paper and plastic. For an adequate comparison paper 
and plastic were tested first as homogeneous materials. Figure 4.11 shows the real and 
imaginary permittivity results for paper, and Figure 4.12 shows the real and imaginary 
permittivity results for plastic. These figures show that paper has a larger permittivity 
and higher loss than plastic. By interspersing layers of paper and plastic an 
inhomogeneous material was constructed. 


Blue Paper: Real Relative Permittivity vs. Frequency 
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Blue Paper: Imaginary Relative Permittivity vs. Frequency 
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Figure 4.11: Real and imaginary results for blue paper 
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Yellow Plastic: Real Relative Permittivity vs. Frequency 



Yellow Plastic: Imaginary Relative Permittivity vs. Frequency 
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Figure 4.12: Real and imaginary results for yellow plastic 
The inhomogeneous material was fashioned using five sub-layers, each with an 
increasingly larger number of paper layers until the last sub-layer contained only paper. 
Thus a material with an initial permittivity intercept and an increasing permittivity profile 
was created. Figure 4.13 shows the construction of the inhomogeneous material. 

The material was fashioned so that it would fit securely inside the sample holder 
in the x and y directions. It was also clamped overnight to ensure that the layers were 
sandwiched tightly and the material’s thickness in the z direction nearly uniform. 

The material was tested so that the yellow plastic portion of the material began at 
the z=0 interface. And the material progressively increased in permittivity in the positive 
z direction (contains more and more blue paper layers). Figure 4.14 shows the results for 
the z-intercept. The extracted z-intercept from the inhomogeneous algorithm 



Figure 4.13: Plastic and paper inhomogeneous material construction 
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is shown in red and is compared to the permittivity of yellow plastic when calculated 
alone using the homogenous algorithm in blue. Both the real and imaginary parts of the 
relative intercept calculated by the inhomogeneous algorithm are fairly accurate. 
However, measurement error does not account for the discrepancy in the imaginary part 
of the permittivity intercept. One also sees that error plays a smaller role in the 
continuous algorithm when compared to the discrete algorithm. 

Figure 4.15 shows the relative slope calculated by the inhomogeneous algorithm 
in red compared to the rise over run slope calculated from the homogeneous algorithms. 
Again, a strong consistency is seen between the inhomogeneous algorithm and the actual 
inhomogeneous permittivity of the material for both real and imaginary parts. It is also 
seen that the algorithm will converge with a material with a relatively small slope. The 
measurement errors again do not account for the slight divergence from the actual profile. 
These errors could be cause from inconsistencies during the fabrication process, since 
each thin slice of plastic or paper is not perfectly uniform. 


Real Relative Intercept vs. Frequency 


Imaginary Relative Intercept vs. Frequency 
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Figure 4.14: Continuous algorithm relative intercept results 
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Real Relative Slope vs. Frequency 


Imaginary Relative Slope vs. Frequency 
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Figure 4.15: Continuous algorithm relative slope results 
4.4 Discrete and Continuous Algorithm Comparison 

Figures (4.16) and (4.17) show the continuous and discrete algorithm profiles 
compared to the actual profile described in Figure 4.13. In both graphs, the 
inhomogeneous algorithm profile is shown in red, and the actual profile in blue. The left 
graph in Figure 4.16 shows the real part of e r (z) versus z in mm, and the right graph 
shows the imaginary part of e r (z) versus z in mm. All plots are done using the mid-band 
frequency of 10 GHz. The results of the continuous algorithm on the real graph show 
strong agreement at z=0, but the algorithm starts to diverge from the actual profile as 
z—>d(i. The results in the imaginary graph are most accurate in the center of the material, 
but the algorithm diverges from the actual profile as z->0 or z-><7 0 . 

Figure 4.17 shows the accuracy of the discrete algorithm. For this material only 
left justified measurements were taken since the material was too thick to also do right 
justified measurements. The discrete algorithm closely follows the actual profile, except 
that at the edges of each step it becomes less accurate. The large uncertainty possibilities 
are also seen in this figure, especially for the £2 section. 
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Continuous Algorithm: Real Relative Permittivity vs. z 



Continuous Algorithm: Imaginary Relative Permittivity vs. z 



Figure 4.16: Continuous algorithm profile results 


Discrete Algorithm: Real Relative Permittivity vs. z 


Discrete Algorithm: Imaginary Relative Permittivity vs. z 




Figure 4.17: Discrete algorithm profile results 
If one were to use a homogeneous algorithm to find the profile of this same 
inhomogeneous material the results would be similar to what is shown in Figure 4.18 and 
Figure 4.19. These show the permittivity as a function of frequency. In Figure 4.20 and 
Figure 4.21 the permittivity is plotted versus position z within the material at the mid¬ 
band frequency of 10 GHz. 

Despite some of the inaccuracies associated with the discrete and continuous 
profiles, they still better articulate the changing profile within the material. 


57 

































































Homogeneous Algorithm w/ Inhomogeneous Media: Refc^ vs Frequency Homogeneous Algorithm w/ Inhomogeneous Media: Im^ vs Frequency 




Figure 4.18: Inhomogeneous material characterized with homogeneous algorithm 
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Homogeneous Algorithm w/ Inhomogeneous Media: Imfc^ vs z 
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Figure 4.18: Inhomogeneous material profiled with homogeneous algorithm 
The discrete and continuous algorithms have their own benefits and drawbacks. 
The discrete algorithm’s advantages include its ability to characterize a material with a 
simpler theoretical foundation since it relies on homogeneous wave equations. It also 
becomes advantageous if one were to characterize a material that is not linear, since it 
would likely provide a more accurate result in comparison to the continuous profile 
because it splits up the material into sections. This approach also lends itself easily to 
characterizing layered stack material configurations such as a Jaumann absorber. 
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A major disadvantage of the discrete algorithm is the resulting profile averages 
across sections within the material. This produces inaccuracies at each permittivity 
section’s right and left boundaries. Also, the discrete algorithm is more error prone as 
can be seen by the error bars. 

One could potentially change the discrete code so that it produces a finer 
characterization with additional sections, but it was discovered during testing that as the 
section numbers were increased the accuracy of the algorithm began to sharply diminish. 
This could be caused by an inability of the network analyzer to distinguish between the 
phase differences from each shim. A more precise measuring device or additional data 
points would allow for a greater set of independent equations for the Levenberg- 
Marquardt root search. It is also possible that an alternative root search method other than 
Levenberg-Marquardt could provide better results. The convergence of the algorithm 
becomes increasingly temperamental with additional root search dimensions. 

The continuous approach is advantageous over the discrete approach in some 
ways including the fact that it produces a smoother profile characterization. It also 
requires only one measurement, so that it takes less effort and time to collect the 
algorithm’s necessary data points. Also, with one measurement there is a lower 
likelihood of error creeping into the resulting profile. 

A disadvantage of this approach in comparison to the discrete method is that it 
requires a more complex theoretical foundation as seen in Chapter 3. Additionally, this 
approach would do a poorer job than the discrete approach when characterizing a non¬ 
linear inhomogeneous profile. Another disadvantage of the continuous algorithm is that 
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the theoretical development produced equations with an instability as m E ->0 as seen from 
(3.19) or (3.20). When characterizing materials with very low slopes a check is 
necessary within the algorithm to ensure that the equations do not become too large and 
lead to an erroneous profile characterization. 

The continuous approach will not scale well with thick materials since it was 
assumed from the beginning of the development that the material was thin. To 
characterize thicker material’s more complex theoretical profiles need to be assumed 
instead of the linear profile used in (3.15). 

4.5 Summary 

This chapter discussed the equipment used for data collection and the results from 
the discrete and continuous algorithms using several testing scenarios. Both the discrete 
and continuous algorithms produced adequate profile characterizations of the test 
materials. However, in some regions within the material they produced slight 
inaccuracies in comparison to the actual pennittivity profile. 

The discrete and continuous algorithms have advantages and disadvantages 
depending on how they are implemented and what materials are being characterized. 
Some basic understanding of the material being investigated is beneficial in selecting the 
best method to use. In the next chapter the entire study will be summarized. 
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V. Conclusion 


This document discussed the theoretical formulation for developing a discrete and 
continuous algorithm to find inhomogeneous permittivity profiles. It revealed the 
algorithms’ results when testing on various materials and the pros and cons of each. 

The discrete approach offered a theoretically simpler framework since it relied on 
homogeneous wave equations in a piece-wise constant permittivity approach. However, 
there was a significant downside to this approach since it required numerous 
measurements to implement. This can lead to gross error in the results if measurements 
are not taken accurately. 

The continuous approach offered a way of creating a smooth permittivity profile 
using inhomogeneous wave equations. This required only one measurement, but relies 
on a more complex development and a thin material assumption. Unfortunately, this 
continuous algorithm may give misleading results for thicker materials. 

Despite some of the drawbacks associated with either technique, they still offer a 
reasonable solution to inhomogeneous permittivity characterization. Depending on the 
material, one may choose one method or the other. 

5.1 Future Research 

Improvements can be made to each approach presented in this study. For the 
discrete case one could improve the precision of the permittivity profile solutions. This 
could be done by using additional measurement data, a one port calibration technique to 
reduce noise, or an alternative root search method. During the shim calculations a ten 
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degrees phase was used for each shim thickness. If it is possible to use a smaller phase 
spacing with a device capable of distinguishing these phase differences, more 
independent measurements can be included in the algorithm for a more detailed profile 
result. 

The continuous case can also be improved. In this study the material was 
assumed to be thin, and the profile approximated as a linear function. Work could be 
done to investigate permittivity profiles of thicker materials, or ones that have parabolic, 
sinusoidal, or more complicated pennittivity profiles. 

Additionally, both approaches found only inhomogeneous permittivity profiles. 
These developments could be expanded to include penneability inhomogeneity as well. 
These approaches could also be combined in the future to produce an algorithm that 
would solve for a piece-wise linear profile. This would allow the characterization of 
more complicated profiles. 
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Appendix A: Constitutive parameter extraction main code example 


%Written by: Capt Marcus A Sitterly 
%Date: 19 July 11 

%Sample Code for constitutive parameter extraction for one material 

%CONSTANTS USED IN CALCULATIONS 
epsilonnot = (10 A -9)/ (36*pi); 
munot = 4*pi*10 A -7; 
a = 2.286*10 A -2; 
lengthl = 6.64*10 A -3; 

Plexiglas 

%DATA MANIPULATION BEFORE INPUT TO LEVENBERG-MARQUAD ALGORITHM 

%Read in laboratory data from text file with custom function 
[ frequency Sllexp S21exp S12exp S22exp ] = 
dataextract( 'plexi thick.txt'); 

omega = 2*pi*frequency; %Angular Frequency 

kO = omega*sqrt(epsilonnot*munot); %Propogation constant k 
kzO = sqrt((kO. A 2)-(pi/a) A 2); %Propogation constant kz 


%F/m 

%H/m 

%'a' dimension in m for WR90 
%length of sample of thick 


%Measurements derived from lab with phase shift 
Sllexp = Sllexp; 

S12exp = S12exp.*exp(-j.*kzO.*lengthl); 

S21exp = S21exp.*exp(-j.*kzO.*lengthl); 

S22exp = S22exp.*exp(-2.*j.*kzO.*lengthl); 

%Allocate space for epsilonrell and murell vectors 
epsilonrell = zeros(length(frequency),1); 
murell = zeros(length(frequency),1); 

%CALCULATE EPSILON RELATIVE AND MU RELATIVE FOR GIVEN FREQUENCY 
for t=l:length(frequency) 

%Experimental values derived from lab [Sll S12 S21 S22] for Plexiglas 
y =[Sllexp(t) S12exp(t) S21exp(t) S22exp(t)]; 

%Vector that contains frequency and length measurement for algorithm 
x=[omega(t) lengthl]; 

%Initial guess [epsilonrell murell] 
betanot= [3,1]; 

%Use NLFIT to obtain epsilon/mu 
beta=nlinfit(x,y,Sonematerial,betanot); 
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%Assign epsilon and mu values from numerical calculation 
epsilonrell(t)=beta(1); 
murell(t)=beta (2); 

end 

%PLOT RESULTS 

%Plot real(e) vs frequency 
figure(1) ; 

plot ( (10 A -9)*frequency, (real(epsilonrell))); 

axis([min((10 A -9)*frequency) max((10 A — 9)*frequency) 0 3]); % axis([xmin 
xmax ymin ymax]) 

title (' FREQUENCY VS. REAL \epsilon_r', 'fontsize', 12) 
xlabel (' Frequency (GHz)', 'fontsize', 12); 
ylabel (' Re(\epsilon_r) ' , 'fontsize', 12); 

%Plot imag(e) vs frequency 
figure (2); 

plot((10 A — 9)*frequency, (imag(epsilonrell))) ; 

axis([min((10 A -9)*frequency) max((10 A — 9)*frequency) -1 1]); % 
axis([xmin xmax ymin ymax]) 

title (' FREQUENCY VS. IMAGINARY \epsilon_r', 'fontsize', 12) 
xlabel (' Frequency (GHz)', 'fontsize', 12); 
ylabel (' Im(\epsilon_r) ' , 'fontsize', 12); 

%Plot real(mu) vs frequency 
figure(3); 

plot((10 A — 9)*frequency, (real(murell))); 

axis([min((10 A -9)*frequency) max((10 A -9)*frequency) 02]); % axis([xmin 
xmax ymin ymax]) 

title (' FREQUENCY VS. REAL \mu_r' , 'fontsize', 12) 
xlabel (' Frequency (GHz)', 'fontsize', 12); 
ylabel (' Re(\mu_r)' , 'fontsize', 12); 

%Plot imag(mu) vs frequency 
figure(4); 

plot ( (10 A — 9)*frequency, (imag(murell))); 

axis([min((10 A -9)*frequency) max((10 A — 9)*frequency) -1 1]); % 
axis([xmin xmax ymin ymax]) 

title (' FREQUENCY VS. IMAGINARY \mu_r' , 'fontsize', 12) 
xlabel (' Frequency (GHz)', 'fontsize', 12); 
ylabel (' Im(\mu_r)' , 'fontsize', 12); 
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Appendix B: Code for nlinfit function 


function y=onematerial(betanot,x) 


a = 2.286*10 A -2; 
omega = x(1); 

epsilonnot = (10 A -9)/(36*pi); 
munot = 4*pi*10 A -7; 


%'a' dimension in m for WR-90 

%Angular Frequency 

%F/m 

%H/m 


o 


MATERIAL SLICES PARAMETERS 


%Material 

epsilonrell = betanot(l); 
epsilon 

murell = betanot(2); 


%betanot(l) is algorithms value for 


%betanot(2) is algorithms value for mu 


lengthl = x(2); 

%Last interface (free space) 
epsilonrel2 = 1; 
murel2 = 1; 
lengthend = 0; 


% PROROGATION CONSTANT CALCULATIONS 

%Propagation constant k 

kO = omega*sqrt(epsilonnot*munot); 

kl = omega*sqrt(epsilonrell*epsilonnot*murell*munot); 
k2 = omega*sqrt(epsilonrel2*epsilonnot*murel2*munot); 

%Propagation constant kz 
kzO = sqrt((kO A 2)-(pi/a) A 2); 
kzl = sqrt((kl A 2)-(pi/a) A 2); 
kz2 = sqrt((k2 A 2)-(pi/a) A 2); 

%IMPEDANCE CALCULATIONS 
%Material impedance calculations 
etaO = sqrt(munot/epsilonnot) ; 

etal = sqrt((murell*munot)/(epsilonrell*epsilonnot)); 
eta2 = sqrt((murel2*munot)/(epsilonrel2*epsilonnot)) ; 

%Wave impedance calculations 
waveimpedenceO = (eta0*k0)/kzO; 
waveimpedencel = (etal*kl)/kzl; 
waveimpedence2 = (eta2*k2)/kz2; 

R1 = (waveimpedencel - waveimpedenceO)/(waveimpedencel + 
waveimpedenceO); 

R2 = (waveimpedence2 - waveimpedencel)/(waveimpedence2 + 
waveimpedencel); 
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PI = exp(-j*kzl*lengthl); 

P2 = exp(-j*kz2*lengthend); 

%A MATRIX CALCULATIONS 

A1 = (1/(PI*(1+R1)))*[1 (R1*P1 A 2) ;R1 P1 A 2]; 

A2 = (1/(P2*(1+R2)))*[1 (R2*P2 A 2);R2 P2 A 2]; 


%Find A Matrix for the system 
ASys = A1*A2; 

All = ASys (1,1); 

A12 = ASys (1,2); 

A21 = ASys(2,1); 

A22 = ASys(2,2); 

%Find the S matrix from the A matrix 

SSys = (1/All)*[A21 (A11*A22-A21*A12);1 -A12] 

511 = SSys(1,1); 

512 = SSys(1,2); 

521 = SSys(2,1); 

522 = SSys(2,2); 

%Compare to labratory data 
y = [Sll S12 S21 S22] ; 
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